---
title: "7-5 Days alive and free of invasive or non-invasive ventilation"
description: |
This outcome was defined as the number of days alive and free of
invasive or non-invasive ventilation from randomisation to 28 days.
All patients dying within 28 days will be assigned zero free days.
author:
- name: James Totterdell
affiliation: University of Sydney
- name: Rob Mahar
affiliation: University of Melbourne
date: today
toc-depth: 5
freeze: true
---
# Preamble
```{r}
#| label: pkgs
#| code-summary: Load packages
library (tidyverse)
library (labelled)
library (kableExtra)
library (cmdstanr)
library (posterior)
library (bayestestR)
library (bayesplot)
library (matrixStats)
library (plotly)
library (lubridate)
library (ggdist)
library (patchwork)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: data
#| code-summary: Prepare dataset
devtools:: load_all ()
all_dat <- read_all_no_daily ()
all_dat_daily <- read_all_daily ()
acs_itt_dat <- all_dat %>%
filter_acs_itt () %>%
transmute_model_cols_grp_aus_nz ()
acs_itt_nona_dat <- acs_itt_dat %>%
filter (! is.na (out_dafv))
# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
filter_concurrent_intermediate ()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>%
filter (! is.na (out_dafv))
# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
filter_concurrent_std_aspirin ()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>%
filter (! is.na (out_dafv))
# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
filter_concurrent_therapeutic ()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>%
filter (! is.na (out_dafv))
```
```{r}
#| label: models
#| code-summary: Load models
ordmod0 <- cmdstan_model (
"stan/ordinal/logistic_cumulative.stan" ) # No epoch or site
ordmod_epoch <- cmdstan_model (
"stan/ordinal/logistic_cumulative_epoch.stan" ) # Epoch only
ordmod <- cmdstan_model (
"stan/ordinal/logistic_cumulative_epoch_site.stan" ) # Full model
```
```{r}
make_summary_table <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" , "C3" , "C4" ),
labels = intervention_labels2 ()$ CAssignment[- 1 ])) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafv)),
Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` Any ventilation ` = sprintf ("%i (%.0f%%)" ,
sum (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE ),
100 * mean (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE )),
` DAFV, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafv, na.rm = T),
quantile (out_dafv, 0.25 , na.rm = TRUE ),
quantile (out_dafv, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafv)),
Deaths = sprintf ("%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` Any ventilation ` = sprintf ("%i (%.0f%%)" ,
sum (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE ),
100 * mean (D28_OutcomeDaysFreeOfVentilation < 28 , na.rm = TRUE )),
` DAFV, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafv, na.rm = T),
quantile (out_dafv, 0.25 , na.rm = TRUE ),
quantile (out_dafv, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Anticoagulation \n intervention ` = CAssignment)
kable (
tdat,
format = format,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
make_primary_model_data <- function (
dat,
vars = NULL ,
beta_sd_var = NULL ,
ctr = contr.orthonorm,
p_mult = 2 ) {
X <- make_X_design (dat, vars, ctr)
attX <- attr (X, "contrasts" )
X <- X[, - 1 ]
attr (X, "contrasts" ) <- attX
nXtrt <- sum (grepl ("rand" , colnames (X)))
beta_sd <- c (rep (1 , nXtrt), beta_sd_var)
epoch <- dat$ epoch
M_epoch <- max (dat$ epoch)
region <- dat$ ctry_num
M_region <- max (region)
site <- dat$ site_num
M_site <- max (site)
region_by_site <- dat %>%
dplyr:: count (ctry_num, site_num) %>%
pull (ctry_num)
y_raw <- ordered (dat$ out_dafv)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
list (N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
M_region = M_region,
M_site = M_site, site = site,
M_epoch = M_epoch, epoch = epoch,
region_by_site = region_by_site,
p_par = p_mult * rep (1 / unique_y, unique_y),
beta_sd = beta_sd)
}
```
## Outcome Derivation
This outcome was defined as the number of days alive and free of positive pressure ventilation
(invasive or non-invasive) to 28 days. All patients dying within 28 days will be assigned zero free days.
The outcome is calculated for a patient as:
- missing, if day 28 mortality is unknown (`D28_PatientStatus` ) or if the patient was known to have been alive
at day 28 but number of days requiring ventilation (`D28_OutcomeDaysFreeOfVentilation` ) is unknown
- 0 if they died by day 28 (`D28_PatientStatus` )
- `min(28, D28_OutcomeDaysFreeOfVentilation)` otherwise
## Descriptive Analyses
### Overall
The overall distribution of days alive and free of ventilation (DAFV) are reported in @fig-dafv-dist for all participants in the ACS-ITT set.
There were few participants who required any ventilation during hospital who did not die (DAFV of 0),
```{r}
#| label: fig-dafv-dist
#| fig-cap: |
#| Distribution of days alive and free of hospital.
#| fig-subcap: ACS-ITT
#| fig-height: 4
pdat <- acs_itt_dat %>%
dplyr:: count (dafv = addNA (ordered (out_dafv, levels = 0 : 28 )), .drop = F) %>%
mutate (p = n / sum (n))
ggplot (pdat, aes (dafv, p)) +
geom_col () +
labs (
x = "Days alive and free of ventilation" ,
y = "Proportion"
)
```
```{r}
#| tbl-cap: Summaty of DAFV.
save_tex_table (
make_summary_table (acs_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-5-anticoagulation-summary" ))
make_summary_table (acs_itt_dat)
```
### Anti-coagulation
```{r}
#| label: fig-dafh-dist-anticoagulation
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment[- 1 ], "<br>" , " \n " )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive.pdf" ), p, height = 3 , width = 6 )
p
```
### Age
```{r}
#| label: fig-dafv-dist-age
#| code-summary: DAFV by age
#| fig-cap: |
#| Distribution of days alive and free of ventilation by age group
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
agegrp = cut (AgeAtEntry, c (18 , seq (25 , 75 , 5 ), 100 ), include.lowest = T, right = F),
dafv = fct_rev (ordered (out_dafv, levels = 0 : 28 )),
.drop = F) %>%
group_by (agegrp) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (agegrp) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (agegrp, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" ) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 ))
p2 <- ggplot (pdat, aes (agegrp, p)) +
geom_col (aes (fill = dafv)) +
labs (x = "Age" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-age.pdf" ), p, height = 2.5 , width = 6 )
p
```
### Country
```{r}
#| label: fig-dafv-dist-country
#| code-summary: DAFV by country
#| fig-cap: |
#| Distribution of days alive and free of ventilation by country.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
country,
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (country) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (country) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (pdat, aes (country, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-country.pdf" ), p, height = 2.5 , width = 6 )
p
```
### Site
```{r}
#| label: fig-dafv-dist-site
#| code-summary: DAFV by site
#| fig-cap: |
#| Distribution of days alive and free of ventilation by study site
#| fig-height: 4
pdat <- all_dat %>%
filter_acs_itt () %>%
filter (! is.na (out_dafv)) %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
Site = fct_infreq (Location),
dafv = ordered (out_dafv, levels = 0 : 28 )) %>%
complete (dafv = ordered (0 : 28 ), nesting (Country, Site), fill = list (n = 0 )) %>%
group_by (Country, Site) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup () %>%
mutate (
Country = droplevels (Country),
Site = droplevels (Site)
)
pdat2 <- pdat %>%
group_by (Country, Site) %>%
summarise (n = sum (n)) %>%
ungroup ()
p1 <- ggplot (pdat2, aes (Site, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (pdat, aes (Site, p)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 1 )) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 / p2 +
plot_layout (guides = 'collect' )
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
### Calendar Time
```{r}
#| label: fig-dafv-dist-calendar
#| code-summary: DAFV by calendar date
#| fig-cap: |
#| Distribution of days alive and free of ventilation by calendar time.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
yr = year (RandDate), mth = month (RandDate),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (yr, mth) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p1 <- pdat %>%
group_by (yr, mth) %>%
summarise (n = sum (n)) %>%
ggplot (., aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p2 <- ggplot (pdat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Calendar date (month of year)" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFV" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" )) +
scale_x_continuous (breaks = 1 : 12 )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-5-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
### Sample Cumulative Logits
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits.
trt_counts <- acs_itt_nona_dat %>%
dplyr:: count (CAssignment, out_dafv) %>%
complete (CAssignment, out_dafv, fill = list (n = 0 )) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (CAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_dafv) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (out_dafv != 28 )
ggplot (trt_logit, aes (CAssignment, rel_clogit)) +
facet_wrap ( ~ out_dafv) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
ggplot (trt_logit, aes (out_dafv, rel_clogit)) +
facet_wrap ( ~ CAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
## Modelling
The primary analysis uses the ACS-ITT set, so all participants have been randomised to the anticoagulation domain and some participants may have been randomised to the antiviral domain.
### Primary Model
```{r}
#| label: fit-model-primary
#| code-summary: Fit primary model
fit_primary_model <- function (dat, ctr = contr.orthonorm, save = FALSE ) {
mdat <- make_primary_model_data (
dat, c ("inelgc3" , "agegte60" , "ctry" ),
c (10 , 2.5 , 1 , 1 ),
ctr)
snk <- capture.output (
mfit <- ordmod$ sample (
data = mdat,
seed = 3510392 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0 ,
adapt_delta = 0.98
))
if (save) {
mfit$ save_object (file.path ("outputs" , "models" , "secondary" , "7-5.rds" ))
}
mdrws <- as_draws_rvars (mfit$ draws ())
# Add names
epoch_map <- dat %>% dplyr:: count (epoch, epoch_lab)
site_map <- dat %>% dplyr:: count (site_num, site)
names (mdrws$ beta) <- colnames (mdat$ X)
names (mdrws$ gamma_epoch) <- epoch_map$ epoch_lab
names (mdrws$ gamma_site) <- site_map$ site
# Convert to treatment log odds ratios
mdrws$ Acon <- attr (mdat$ X, "contrasts" )$ randA %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]
mdrws$ trtA <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]
mdrws$ trtC <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c (
"Intermediate vs low" = exp (mdrws$ trtC[1 ]),
"Low with aspirin vs low" = exp (mdrws$ trtC[2 ]),
"Therapeutic vs low" = exp (mdrws$ trtC[3 ]),
"Intermediate vs low with aspirin" = exp (mdrws$ trtC[1 ] - mdrws$ trtC[2 ]),
"Intermediate vs therapeutic" = exp (mdrws$ trtC[1 ] - mdrws$ trtC[3 ]),
"Low with aspirin vs therapeutic" = exp (mdrws$ trtC[2 ] - mdrws$ trtC[3 ])
)
mdrws$ OR <- c (
setNames (exp (mdrws$ trtC), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),
"Ineligible aspirin" = exp (mdrws$ beta[which (grepl ("inelg" , colnames (mdat$ X)))]),
"Age \u2265 60" = exp (mdrws$ beta[which (grepl ("age" , colnames (mdat$ X)))]),
setNames (exp (mdrws$ beta[which (grepl ("ctry" , colnames (mdat$ X)))]), c ("AU/NZ" , "NP" ))
)
return (list (dat = mdat, fit = mfit, drws = mdrws))
}
res <- fit_primary_model (acs_itt_nona_dat, save = TRUE )
```
```{r}
#| label: odds-ratio-summary-table-primary-model-save
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects).
save_tex_table (
odds_ratio_summary_table_rev (res$ drws$ OR, "latex" ),
"outcomes/secondary/7-5-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table_rev (res$ drws$ OR)
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-5-primary-model-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
```{r}
#| fig-cap: Posterior densities for treatment comparisons.
p <- plot_or_densities (res$ drws$ compare)
p
```
:::{.callout-caution collapse="true"}
###### Diagnostics
```{r}
res$ fit$ summary (c ("alpha" , "beta" , "gamma_site" , "gamma_epoch" , "tau_site" , "tau_epoch" )) %>%
print (n = Inf )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
###### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["alpha" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
#### Posterior Predictive
```{r}
y_raw <- as.integer (levels (res$ dat$ y_raw))
y_ppc <- res$ drws$ y_ppc
y_ppc_raw <- rfun (\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols (acs_itt_nona_dat, tibble (y_ppc = y_ppc_raw)) %>%
mutate (CAssignment = factor (
CAssignment,
levels = names (intervention_labels_short_break ()$ CAssignment),
labels = intervention_labels_short_break ()$ CAssignment))
grp_ppc <- function (grp, d = 0 : 28 ) {
ppc_dat %>%
group_by (grp = {{grp}}) %>%
summarise (
y_lteq = map_dbl (d, ~ mean (out_dafv <= .x)),
ypp_lteq = map (d, ~ rvar_mean (y_ppc <= .x)),
y_eq = map_dbl (d, ~ mean (out_dafv == .x)),
ypp_eq = map (d, ~ rvar_mean (y_ppc == .x))
) %>%
unnest (c (ypp_lteq, ypp_eq)) %>%
mutate (days = d) %>%
ungroup () %>%
pivot_longer (
y_lteq: ypp_eq,
names_to = c ("response" , "event" ),
names_sep = "_" ,
values_to = "posterior" )
}
plot_grp_ppc <- function (dat, lab = "" , xlab = "Probability" ) {
ggplot (dat %>%
filter (response == "ypp" , event == "eq" ),
aes (y = days)) +
facet_wrap ( ~ grp, nrow = 1 , scales = "free_x" ) +
stat_slabinterval (aes (xdist = posterior), fatten_point = 1 ) +
geom_point (data = dat %>% filter (response == "y" , event == "eq" ),
aes (y = days, x = mean (posterior)),
colour = "red" ,
shape = 23 ) +
scale_x_continuous (
xlab, breaks = c (0 , 0.3 , 0.6 ),
sec.axis = sec_axis (
~ . , name = lab, breaks = NULL , labels = NULL )) +
scale_y_continuous ("Days" ) +
theme (strip.text = element_text (size = rel (0.7 )),
axis.title.x = element_text (size = rel (0.7 )),
axis.text.x = element_text (size = rel (0.65 )),
axis.title.y = element_text (size = rel (0.75 )),
axis.title.x.bottom = element_blank ())
}
pp_C <- grp_ppc (CAssignment)
pp_ctry <- grp_ppc (country)
pp_epoch <- grp_ppc (epoch)
pp_site <- grp_ppc (site) %>%
left_join (ppc_dat %>% dplyr:: count (site, country),
by = c ("grp" = "site" ))
p1 <- plot_grp_ppc (pp_C, "Anticoagulation" , "" )
p2 <- plot_grp_ppc (pp_ctry, "Country" , "" )
p3 <- plot_grp_ppc (pp_epoch, "Epoch" , "" )
p4 <- plot_grp_ppc (
pp_site %>% filter (country == "IN" ), "Sites India" , "" )
p5 <- plot_grp_ppc (
pp_site %>% filter (country == "AU" ), "Sites Australia" , "" )
p6 <- plot_grp_ppc (
pp_site %>% filter (country == "NP" ), "Sites Nepal" , "" )
p7 <- plot_grp_ppc (
pp_site %>% filter (country == "NP" ), "Sites New Zealand" , "" )
g1 <- (p1 | p2) / p3
g2 <- p4 / p5 / (p6 | p7)
pth1 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-5-primary-model-acs-itt-ppc1.pdf" )
pth2 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-5-primary-model-acs-itt-ppc2.pdf" )
ggsave (pth1, g1, width = 6 , height = 5 , device = cairo_pdf)
ggsave (pth2, g2, width = 6 , height = 6 , device = cairo_pdf)
g1
g2
```
### Sensitivity: Concurrent controls
Three separate analyses are conducted based on concurrent randomisations as was done for the primary outcome. For these models, the epoch term has been removed.
##### Intermediate dose
```{r}
#| code-summary: Descriptive summary table
save_tex_table (
make_summary_table (acs_itt_concurc2_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-5-anticoagulation-concurrent-intermediate-summary" ))
make_summary_table (acs_itt_concurc2_dat)
```
```{r}
#| code-summary: Fit model
ctr <- contr.orthonorm
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc2_nona_dat,
contrasts.arg = list (CAssignment = ctr))[, - 1 ]
y_raw <- ordered (acs_itt_concurc2_nona_dat$ out_dafv)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
mdat <- list (
N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
p_par = 2 * rep (1 / unique_y, unique_y),
beta_sd = c (1 , 2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- ordmod0$ sample (
data = mdat,
seed = 75136 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mdrws <- as_draws_rvars (mfit$ draws ())
names (mdrws$ beta) <- colnames (X)
mdrws$ Ccon <- as.vector (ctr (2 ) %**% mdrws$ beta[1 ])
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]))
mdrws$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt),
setNames (exp (mdrws$ beta[2 : 4 ]), c ("Age \u2265 60" , "Australia/New Zealand" , "Nepal" )))
```
```{r}
#| code-summary: Sample distribution
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention in subset.
#| fig-height: 4
pdat <- acs_itt_concurc2_nona_dat %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" ),
labels = c ("Low" , "Intermediate" )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv)), colour = NA ) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive-concurrent-intermediate.pdf" ), p, height = 2.5 , width = 4 )
p
```
```{r}
#| code-summary: Posterior contrast
#| tbl-cap: Posterior densities for comparisons of interest.
p <- plot_or_densities (mdrws$ compare) +
geom_text (aes (x = 1 , y = 1 - 0.1 ,
label = sprintf ("Pr(OR > 1) = %.2f" , Pr (RV > 1 ))),
hjust = 0 , nudge_x = 0.01 , family = "Palatino" ) +
scale_x_log10 (breaks = c (0.8 , 0.9 , 1 , 1.1 , 1.2 , 1.3 , 1.4 , 1.5 , 1.6 )) +
labs (y = "Comparison" )
p
```
```{r}
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),
"outcomes/secondary/7-5-primary-model-acs-itt-concurrent-intermediate-summary-table" )
odds_ratio_summary_table_rev (mdrws$ OR)
```
:::{.callout-caution collapse="true"}
###### Diagnostics
```{r}
mfit$ summary ()
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
###### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["alpha" ])
```
:::
#### Low-dose with aspirin
```{r}
#| code-summary: Descriptive summary table
save_tex_table (
make_summary_table (acs_itt_concurc3_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-5-anticoagulation-concurrent-stdaspirin-summary" ))
make_summary_table (acs_itt_concurc3_dat)
```
```{r}
#| code-summary: Fit model
ctr <- contr.orthonorm
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc3_nona_dat,
contrasts.arg = list (CAssignment = ctr))[, - 1 ]
y_raw <- ordered (acs_itt_concurc3_nona_dat$ out_dafv)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
mdat <- list (
N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
p_par = 2 * rep (1 / unique_y, unique_y),
beta_sd = c (1 , 1 , 2.5 , 1 )
)
snk <- capture.output (
mfit <- ordmod0$ sample (
data = mdat,
seed = 135356 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mdrws <- as_draws_rvars (mfit$ draws ())
names (mdrws$ beta) <- colnames (X)
mdrws$ Ccon <- as.vector (ctr (3 ) %**% mdrws$ beta[1 : 2 ])
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),
"Low with aspirin vs low" = exp (mdrws$ Ctrt[2 ]),
"Intermediate vs low with aspirin" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))
mdrws$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt[1 ]),
"Low with aspirin" = exp (mdrws$ Ctrt[2 ]),
setNames (exp (mdrws$ beta[3 : 4 ]), c ("Age \u2265 60" , "Australia/New Zealand" )))
```
```{r}
#| code-summary: Sample distribution
#| fig-cap: |
#| Distribution of days alive and free of ventilation by anti-coagulation intervention in subset.
#| fig-height: 4
pdat <- acs_itt_concurc3_nona_dat %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" , "C3" ),
labels = c ("Low" , "Intermediate" , "Low \n with aspirin" )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) +
theme (legend.key.size = unit (0.75 , "lines" ))
p
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive-concurrent-stdaspirin.pdf" ), p, height = 2.5 , width = 6 )
```
```{r}
#| code-summary: Posterior contrast
#| tbl-cap: Decision summaries on comparisons
plot_or_densities (mdrws$ compare)
```
```{r}
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),
"outcomes/secondary/7-5-primary-model-acs-itt-concurrent-stdaspirin-summary-table" )
odds_ratio_summary_table_rev (mdrws$ OR)
```
:::{.callout-caution collapse="true"}
###### Diagnostics
```{r}
mfit$ summary ()
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
###### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["alpha" ])
```
:::
#### Therapeutic dose
```{r}
#| code-summary: Descriptive summary table
save_tex_table (
make_summary_table (acs_itt_concurc4_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-5-anticoagulation-concurrent-therapeutic-summary" ))
make_summary_table (acs_itt_concurc4_dat)
```
```{r}
#| code-summary: Fit model
ctr <- contr.orthonorm
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc4_nona_dat,
contrasts.arg = list (CAssignment = ctr))[, - 1 ]
y_raw <- ordered (acs_itt_concurc4_nona_dat$ out_dafv)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
mdat <- list (
N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
p_par = 2 * rep (1 / unique_y, unique_y),
beta_sd = c (1 , 1 , 2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- ordmod0$ sample (
data = mdat,
seed = 49135 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mdrws <- as_draws_rvars (mfit$ draws ())
names (mdrws$ beta) <- colnames (X)
mdrws$ Ccon <- as.vector (ctr (3 ) %**% mdrws$ beta[1 : 2 ])
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),
"Therapeutic vs low" = exp (mdrws$ Ctrt[2 ]),
"Intermediate vs therapeutic" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))
mdrws$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt[1 ]),
"Therapeutic" = exp (mdrws$ Ctrt[2 ]),
setNames (exp (mdrws$ beta[3 : 5 ]), c ("Age \u2265 60" , "India" , "Australia/New Zealand" )))
```
```{r}
#| code-summary: Sample distribution
#| fig-cap: |
#| Distribution of days alive and free of ventilation by anti-coagulation intervention in subset.
#| fig-height: 4
pdat <- acs_itt_concurc4_nona_dat %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" , "C4" ),
labels = c ("Low" , "Intermediate" , "Therapeutic" )),
dafv = ordered (out_dafv, levels = 0 : 28 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafv))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of ventilation" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-5-descriptive-concurrent-therapeutic.pdf" ), p, height = 2.5 , width = 6 )
p
```
```{r}
#| code-summary: Posterior contrasts
#| tbl-cap: Decision summaries on comparisons
plot_or_densities (mdrws$ compare)
```
```{r}
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),
"outcomes/secondary/7-5-primary-model-acs-itt-concurrent-therapeutic-summary-table" )
odds_ratio_summary_table_rev (mdrws$ OR)
```
:::{.callout-caution collapse="true"}
###### Diagnostics
```{r}
mfit$ summary ()
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
###### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["alpha" ])
```
:::